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Abstract 



Nonlinear systems of the reaction-diffusion type, including Gierer-Meinhardt 
models of autocatalysis, are studied by using Lie algebras coming from the 
prolongation structure. The consequences of this analytical approach, as the 
determination of special exact solutions, are compared with the corresponding 
results obtained via numerical simulations. 
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Typeset using REVTf^X 



I. INTRODUCTION 



Models for biological pattern formation may be described by reaction-diffusion (RD) 
equations of the type []J 

u k t = a k V 2 u k + R k {u), (1.1) 

(no sum over k), where u = {u k } are the dynamical fields, k = 1,2, . . . , N, a k are constants 
and Rk{u) are functions of the fields which characterize the reactions among them. 

The class of Eq. (\1 . 1| ) includes popular models of chemical reactions |§. In particular, 
when the diffusion coefficients a k are not all equal, Turing instabilities may occur and com- 
plicated patterns emerge which have been related to morphogenesis, reaction front dynamics, 
self-organization and so forth 0. 

Models of pattern generation in complex organisms are investigated mostly via computer 
simulations. A minor attention is payed to the analysis of the mathematical properties of the 
underlying equations. The knowledge of these properties, combined with numerical studies, 
could be of help both to compare theory and experiment, to check the validity of the models, 
and to suggest possible improvements. 

Keeping in mind this programme, in this paper we apply the prolongation technique || 
to the 1 + 1 dimensional RD system 

u kt = a k u kxx + R k (u), (1.2) 

(k = 1,2, ... , N), where first the reaction functions R k (u) are not specified. Our approach 
contemplates the use of the concept of pseudopotential, that is an M-component vector (M 
arbitrary) y = {y(x, t)} = (y 1} y M ) defined by 

yx = f J (u)T j y, (1.3a) 
y t = g J (u,u x )Tjy, (1.3b) 

where summation over repeated indices is understood, u = {u k }, u x = [u kx ], Tj are M x M 
matrices spanning a Lie algebra £ 

[T l ,T 3 } = c k l3 T k , (1.4) 

and f\ gi are functions to be determined in such a way that the compatibility condition 
Utx = Vxt reproduces Eq. (|L2|). When this feature occurs, then Eqs. ( |1.3| ) represent a 
linearization of Eq. ( |1.2| ). 

The motivation for dealing with Eq. ( |1.2|) instead of Eq. (|1 . 1|) within the prolongation 
scheme is at least threefold. First, the prolongation machinery in 2 + 1 dimensions is not yet 
well established, even if some examples exist concerning a few integrable cases ||. Second, 
biological structures modeled by equations of the class ( |1.2| ) are interesting also in 1 + 1 
dimensions. One of them, which consists of Eqs. (|1.2|) where R k (u) is assumed to be linear 
in the fields {u k } (k = 1,2), has been discussed recently by Kondo and Asai in relation to 
the stripe pattern mechanism of the angelfish Pomachantus 0. Third, the boundary of a 
complex two-dimensional structure hosts a restricted 1 + 1 dimensional RD system || . 

The main results achieved in this work are listed in correspondence of the contents of 
the different Sections in which the paper is planned. 
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In Sec. H, we deal with the prolongation of RD systems described by Eq. ( \1.2\) . In Sec. |T1 



we study a class of quadratic models which are completely linearizable and provide several 
exact analytical solutions of the travelling wave type. These solutions are obtained by direct 
inspection. In Sec. |TV| , we emphasize the role of linearizability. The evolution equations 
for the pseudopotential are written explicitly and solved showing how new solutions can 
emerge in terms of known ones. In Sec. |V], we apply the general results found in Sec. |II 
to the prolongation of the Gierer-Meinhardt models M. The resulting algebra turns out 



to be closed and is that of the similitude group in the plane M. In Sec. VI, we explore 



the consequences of this underlying algebraic structure. In particular, we treat analytical 
approximations of particular classes of solutions, namely homogeneous and of the travelling 
wave type. Section |VII| contains results obtained via numerical integration of the evolution 
equations in order to check the analytical predictions and to suggest possible developments 
of the analysis. Finally, in Sec. |VIII| some concluding remarks are considered, while in the 
Appendices details of the calculations are reported. 



II. PROLONGATION OF GENERAL RD SYSTEMS 

Following the general strategy outlined in the Introduction, we begin our analysis of the 
constraints on the algebra C by writing explicitly the compatibility condition y tx = y x t of 
Eqs. ( |1.3| ), namely 

\ dp dg 1 dg 1 1 

^ { t. — Ukt — t. — u kx — — u kxx > + cZf,f a qP = 0. (2.1) 

V\0Uk du k 8u kx kxx ) a(3J y K J 



Substitution from Eq. ( |1.2|) into Eq. ( [2.1|) yields 

dp . . dg 7 dg 



(/<< i/y^ 

. , {akU kxx + R k ) - ^-u kx - ^—u kxx + clpPg 13 = 0. (2.2) 
du k du k du kx p 

The coefficient of u kxx must vanish 

dp dg^ . . 

du k du kx 

so that 

g 1 = a k ^u kx + A~<(u), (2.4) 
du k 

where A 7 is a function of integration. Putting this result back into Eq. ( |2.2| ), we find 

dP d 2 p dA~< df ,y p 

- — R k - a k - — x—u kx ui x - — — u kx + cUj a k - — u kx + c^ f A p = 0. (2.5) 

du k du k dui du k p du k p 

Consequently, the following three conditions 

6 

du k dui 



d 2 P 

otk + oti)- — — = 0, (2.6a) 



J-R k + d r A^O, (2.6b) 
du k 

<9A 7 dfP 
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hold. Equations Q2.6Q entail two different cases distinguished by the condition 

a k + a^0 Vfc, /, (2.7) 



or 



ot k + oil = for some k, I. (2.8) 

Let us call "non-degenerate" and "degenerate" the problems corresponding to Eq. ( |2.7| ) and 
Eq. (|2.8|), respectively. In these cases, we see that / is linear in the fields u and, therefore, 
the structure of / and g may be expressed by 

f a = H%u k + K a , (2.9a) 
g a = a k H k t u kx + A a (u), (2.9b) 

where H k and K a are constants. Inserting these quantities into Eq. ( |2.6c| ) we obtain 



<9A 7 



clp{H?u x + K a )a k El = a k clpH?Hfa + c^K a H p k a k . (2.10) 



du k 

On the other hand, the compatibility condition 

<9 2 A 7 <9 2 A 7 



du k dui duidu k ' 

gives 



(2.11) 



cl p H^H p k a k = clpHZHfa (a* + a^H^H? = 0. (2.12) 
From this relation, we get 

cZpHZHf = (2.13) 
in the non-degenerate case. Now we can exploit this equation to determine the A 7 term, i.e. 

cLK a Hla k A 7 = cl K a H p k a k + D\ (2.14) 



(D 7 being constants of integration), and the final expressions for / and g, which turn out 
to be 

f a = H k *u k + K a , (2.15a) 

g a = a k H«u kx + c a ^Hla k u k + D a . (2.15b) 

At this stage the reaction terms come into play. All the unknown constants must be chosen 
in order to satisfy the final equation 

^ + C y a A^0. (2.16) 
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To recognize the algebraic structure of the problem, let us adopt a notation in which the 
Lie algebra indices are understood. We define 



/ = f a T a , g = g a T a , (2.17a) 

K = K a T a , D = D a T a , (2.17b) 

$ = H k u k = H^u k T a , (2.17c) 

= a k H k u k = a k H«u k T Q , (2.17d) 



and write the linearized problem in the compact form 



/ = $ + (2.18a) 
g = * x +[K 1 *] + D, (2.18b) 



where 



[H k ,Hi] = 0, (2.19a) 
£ # fc i4 + [$ + AT, [AT, *]+£)] = 0. (2.19b) 

Equations ( |2.19| ) must be satisfied whenever the fields u are chosen. Taking the independent 
monomials in the fields, we obtain an incomplete Lie algebra, in the sense that not all 
the commutators are known. By construction, we know that, if the incomplete algebra is 
satisfied, then the compatibility condition for the linearized problem is assured as long as 
the evolution equations hold. 

On the other hand, the explicit form of the compatibility condition y xt = yt x is easily 
found to be 

$ t - * xx + [$ + K, [K, tf] + D] = 0, (2.20) 

which can be written as 

$t = *xx + y E H k R k ( 2 - 21 ) 

k 

by virtue of Eq. ( |2.19b| ). If the elements {H k } of C were independent, then we could project 
Eq. ( 2.21 ) onto its components and recover the full set of evolution equations estabilishing 
the complete linearization of the reaction-diffusion system under consideration, i.e. 

Vtx = Uxt evolution equations (1.2). (2.22) 

However, as we said before, {H k } are subject to the incomplete Lie algebra (|2.19j ) which is 
a severe constraint. Actually, Eqs. ( |2.19| ) say that every nonlinear term higher than a cubic 
polynomial in Eq. ( |1.2[ ), implies a linear dependence among the elements of the set {H k }. 
In other words, we have 

evolution equations — > y tx = y x t, (2.23a) 
Vtx = Uxt —> linear combinations of evolution equations, (2.23b) 
at most quadratic in the fields. (2.23c) 
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Therefore, the full equivalence expressed by Eq. ( |2.22|) may be achieved only if (i) the reaction 
terms are at most quadratic or (ii) in the special degenerate cases where / and g are no 
more constrained to be linear in the fields. In all the other cases, Eqs. ( |1.3| ) represent only a 
semi-linearization of the evolution equations, i.e. they are equivalent only to a reduced set 
of linear combinations of the original equations. For simplicity, also in these cases we shall 
keep calling Eqs. (|1.3| ) the linearized problem for Eq. (|1.2j) . 

We are interested in applications of the react ion- diffusion systems in problems related to 
morphogenesis, chemical autocatalysis and biological modeling. In these contexts negative 
diffusion coefficients do not admit a simple interpretation and therefore we shall not consider 
in detail the degenerate case. To this regard, we limit ourselves to mention the system 

u it — Uixx + 2u\u2 — 2aui = 0, (2.24a) 
U2t + u 2 xx - 1u\u\ + 2au 2 = 0, (2.24b) 



where a is a fixed constant. Equations ( |2.24j ), which admit an infinite dimensional prolonga- 



tion Lie algebra endowed with a loop structure |J, are integrable, and emerge in the gauge 
formulation of the 1 + 1 dimensional gravity, where u\ and U2 have the meaning of Zweibein 



fields [ 10 1 . These equations are similar to the "fictitious" or "mirror- image" systems with 



negative friction, which appear into the thermo-field approach to the damped oscillator HIT 
Anyway, the role (if any) of Eqs. (|2.24 ) in the modeling of complex organisms remains to 
be elucidated. 



III. LINE ARIZ ABLE QUADRATIC MODELS 

Here we shall build up RD systems which are quadratic and linearizable. We remind the 
reader that some of these models find applications in the study of isothermal autocatalytic 
chemical reactions [13]. In this context it is important the existence of propagating fronts 
(or travelling waves) describing the advance of the reaction. Thus, the investigation of RD 
systems allowing exact solutions of this kind can be a guide for the construction of models 
which may interprete realistic chemical situations. 

Let us start from a closed prolongation algebra and study the explicit form of the com- 
patibility condition expressed by y x t = y x t- In the non-degenerate case, this condition is 
Eq. (|2.20|) , namely a set of quadratic evolution equations, one for each independent Lie al- 
gebra generator appearing after the expansion of the commutators. Since $ and ^ are the 
field dependent terms, the quadratic, linear and constant terms are respectively 

[$, [K, n , [$, D] + [K, [K, n , [K, D] . (3.1) 

Now, let us choose a definite Lie algebra £. First, we fix the set of commuting elements 
H = {H k } involved in the definition of $ and \l/ (see Eq. ( |2. 17(f )). Second, we choose the 
elements K and D. Finally, we evaluate the quantities ( |3.1| ). The contributions proportional 
to the generators H give genuine evolution equations. All the other possible contributions 
fix constraints on the fields. 

If we look for non trivial systems, free of constraints, then we must solve the following 
algebraic problem: find a Lie algebra £ such that: (i) £ has an abelian subalgebra A with 
dimension greater than 2 and (ii) given a basis {Hi} of A there exists K e £ such that 
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[Hi, [Hj, K]] E A. (3.2) 

In theory, we do not know if this problem admits solutions. For instance, if we assume £ 
to be a semisimple algebra and identify A with its Cartan subalgebra, then in the Cartan- 
Weyl basis [n| we have 

[H u Hj] = 0, (3.3a) 
[H u E a ] = aiE a , (3.3b) 

[E a , Ep] = NafiEa+p, (3.3c) 

[E a ,E. a ] = a i H l , (3.3d) 



and for a general K 



we obtain 



^ = E^ + E7«^. (3-4) 



[H it [Hj, K]) = oiiOij~i a E a , (3.5) 

a 

which does not belong to A. 

However, a simple class of algebras with the desidered property is provided by the Ansatz 

[K, D] = 0, [K, A] ~ D, [K, D] e A. (3.6) 

Indeed, application of the Jacobi identity determines all the commutators which turn out to 
be 

[Hi,Hj] =0, i,j = l,...N (3.7a) 

[D,K] = 7i Hi, (3.7b) 

[K, Hi] =^D, (3.7c) 

[D, H^ =fi i X J H j , (3.7d) 

where 7, [i and A are N components vectors which must satisfy 

fi T \ = 0. (3.8) 

Now, let u stand for the column vector of the fields and 

A — diag(ai, - • • , cKjv)- (3.9) 

Expanding the terms in Eq. (|3.1| ) we get 

u t = Au xx + (A yU T u + 7)(/i T v4w + 1). (3.10) 

In the case iV = 2 it is natural to perform a change of variables and introduce in place of ui 
and u 2 , the new fields 

X = fx T u, Y = fi T Au. (3.11) 
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Exploiting the relation /z T A = it is straightforward to show that the two evolution equations 
for (X, Y) are 

X t = Y xx + f i T 1 (X+l), (3.12a) 
Y t = (at + a 2 )Y xx - a x a 2 X xx + AX X + fi T A / y)(Y + 1). (3.12b) 

It is convenient to make the shift 

y->Y + l, (3.13b) 

in Eqs. (|3.12|) . This yields 

X t = Y xx + l x T 1 Y, (3.14a) 

Y t = («i + a 2 )Y xx - a x a 2 X xx + (3(a x - a 2 )XY, (3.14b) 

where 

/j, T AX = Hi\i(a.i — a 2 ) = /?(«! — a 2 ), (3.15) 

P = /xiAi and Eq. (|3.8|) has been used. Now, let us look for solutions of the travelling wave 
type where both X and Y depend on ^ = x + vt. Furthermore, by taking /i T 7 = and, 
consequently, 7 = (X (£ being a constant factor), the shift in Eq. (|3.13|) becomes simply 

X^X + C, (3.16a) 
Y^Y + 1. (3.16b) 

Moreover, we can integrate the first of Eqs. ( |3.14j ) and obtain 

X = -Y' + c , (3.17) 

v 

where Y' = dY/d^ and Co is a constant of integration. The second equation reads 

V Y' = ( ai + a 2 )Y" - ^y'" + ^( ai - a 2 )YY' + /5c F, (3.18) 

v v 

which entails 

Y" - v (— + —)y'+ —Y - — ^-(ai - a 2 )Y 2 + ci = (3.19) 
\ot\ a 2 J ol\Oi 2 Lol\Ol 2 

for Co = 0, with c\ arbitrary constant. 

When the two diffusion constants (ai,a 2 ) take the special values (a, 0) or (a, —a) one 
of the coefficients in Eq. ( |3. 19| ) vanishes. Below, we shall discuss separately these singular 
cases and the general one. 

Anyhow, independently from the value of (a%, a 2 ), we can write the following expressions 
for the fields u\ and u 2 in terms of X and Y (see Eq. ( p.llj) ) and for the relative evolution 
equations 
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"1 



u 2 



1 

/il(« 2 - «i) 

A 2 1 

Ai yUi(a 2 — Oil) 



(a 2 X-Y), 

(atX-Y), // 2 



Ai 



-Hi, 



u u = oi\U\ xx + Ai-R(«i, w 2 ), 

w 2< = a 2 u 2x x + A 2 -R(mi,M2), 



R{ux,u 2 ) 



I aiMi — — a 2 u 2 J + - 



A. Case I: a\ = a, a 2 = 



Equation ( |3.19|) becomes 



v 8 n 

y = -y - — y 2 , 



which may be integrated to give 



Computing X from Eq. (|3.17|) (with Co = 0), we find 



T = tanh— (^-6) 
2a 



2a 2 n±X 



■(1-T 2 )-C, 



y 



a/xiA; 



-(1 + T)-1, 



with the help of Eq. |U6|). 



We have 



B. Case II: a\ = a, a 2 = —a 



Y" - —Y + ^ + ci = 0, 



which furnishes 



Y' 2 = k 1 + k 2 Y + ^- 2 -^-Y 



3 a 



where k± and fc 2 are arbitrary constants. The change of variable 



6a v 2 
Y = — -<p + 



leads to 



/2 A 3 
^ =Atp -g 2 tp- 93, 

-Pi, V " 
^ 2 "fo 2+ 12^' 



93 



1? 



36a 2 



rA'i 



(3v 2 
72a 3 ' 



216a 6 : 



and therefore 



J'(f) = ~n{-6.*,«.) + ^ 



(3.26) 

(3.27a) 
(3.27b) 

(3.27c) 
(3.28) 



where "P is the Weierstrass function [H] and £o an arbitrary (complex) constant. In the 
particular case k\ = k 2 = 0, we have A = g% — 27 g 2 = which involves elementary functions 
only. Using 



P(z, 12c 2 , -8c 3 ) = c + 



3c 



sinh (zy3c) 



we find the following expressions for the fields X, Y 

3v 2 



X 
Y 



2a 2 /xiAi 
3v 2 



Til - T 



2a/^iA] 



(1 — T 



T = tanh— (£-&)■ 



(3.29) 



(3.30a) 

(3.30b) 
(3.30c) 



C. Case III: ai and «2 arbitrary 

In general, Eq. fl3.19|) can be written as 

Y" = aY' + bY + cY 2 + Cl , 
where the coefficients a, b and c are defined by 



a = v 



„,2 



aia2 



2ttiQ!2 

If we introduce the new dependent variable 



(«i - a 2 ). 



(3.31) 

(3.32a) 
(3.32b) 
(3.32c) 
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Y = Y - Y , (3.33) 

Y » = -l{k a2 + h )> ^ 

Cl = -bY - cY 2 , (3.35) 

Y" = aY' - —a 2 Y + cY 2 . (3.36) 
In Appendix |A|, we show that the above equation affords the exact solution 

Y(0 = exp (^A p(lJEexp(4)+ K 0, h) , (3.37) 



where 



and choose c\ to be 



then Y satisfies 



a V 6 y 5 

(ko and k\ are constants of integration) corresponding to the socalled equianharmonic 



case 



14[ . For fco = 0, we deduce 

3a 3 



A ' = 2i^ (1 - T)(1+T)2 - C - (3 ' 38a) 

T = tanh— , (3.38c) 
10' V ' 



where V(z, 0, 0) = z 2 has been exploited. 



D. Case IV: A particular solution for fi T ^f / 



The condition /i T 7 = has been crucial in the previous calculations because it has 
allowed the integration of Eq. ( |3.14aj ). In such a way, some exact solutions turn out to be 
polynomials in tanh(fc£) for a certain k. In the case /z T 7 ^ 0, solutions of the similar kind 
can be obtained by inserting the Ansatz 



N+l 

X = Y, a n T n , 

n=Q 

T = tanh(A;£), 



N 



Y = J2 b n T n , 



n=0 



(3.39a) 
(3.39b) 



into Eqs. (|3.14j) (N arbitrary). The analysis of the case N = 2, provides the solution 



X 



3 

P2V 

3D 



,2\ 3/2 



(3-Dpi - v 2 )\ 5 Pl — —T — D 5pi — — 



D 



D 



Pili - P2I2 
2/iiAi 



Y = 5 Pl - - (1 - T*) - 1 



P2 



D 



T = tanh ( -W 5pi - — 



(3.40a) 
(3.40b) 
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where 



«i + a 2 = 0, D = ai« 2 , pi = /i 7, p 2 = /i AA. (3-41) 

In Fig. [I], we plot the fields Ui(x, 0) and w 2 (x, 0) in the four cases and with a choice of 
the parameters (specified in the captions) such that reaction fronts arise. 

IV. PSEUDOPOTENTIAL FORMULATION OF A QUADRATIC MODEL 

In this Section, we formulate and solve the pseudopotential equations for a particular 
linearizable quadratic model. We show that a bootstrap structure emerges and new solutions 
can be obtained in terms of the old ones found by direct inspection. In doing so, let us assume 

A = (1,1), 7 = (0,0), /i=(l,-l), (4.1) 

in Eq. ( |3.1U| ), which becomes the pair of evolution equations 

u lt = aiu lxx + R(u 1 ,u 2 ), (4.2a) 

u 2t = a 2 u 2xx + R{ui,u 2 ), (4.2b) 

with 

R(ui, u 2 ) = (ui - u 2 )(a 1 u 1 - a 2 u 2 + 1). (4.3) 

In this case, the algebra ( |3~7| ) reads 

[H U H 2 ] = [D, K] = 0, (4.4a) 
[K,^] =-[K,H 2 ] = D, (4.4b) 
[D, Hi] = - [D, H 2 ] = H 1 + H 2 . (4.4c) 

Up to now we have always interpreted the abstract elements Hi, H 2 , K and D as matrices 
belonging to a given ^-dimensional linear representation of the algebra. In such a case, the 
evolution equations for the pseudopotential (see Eqs. ( |2.17|) and (|2.18|) ) take the form 



Vix = Fij(u)yj, (4.5a) 

Vix = Giji^u^yj, (4.5b) 

where F^- and Cry are field dependent matrices. Hereafter, for convenience, we shall write 
Eqs. (|4.5|) in the operator form 



y x = T, (4.6a) 

y t = q, (4.6b) 

where 

y = yidi, T = Fiji^yjdi, G = Giji^yjdi, (4.7) 
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with di = d/dyi. In this way, one has associated with each abstract element Hi, H 2 , K 
and D a differential operator (vector field) whose components, in the basis {di}, are linear 
functions of the pseudopotential variables. It is easy to see that this limitation is not 
necessary. Actually, all the equations go unchanged if arbitrary vector fields are considered. 
Taking this more general attitude, from the relation 

[T, Q\=TQ-QT= -\F, GlijVjdi, (4.8) 

we see that the vector fields Hi, H 2 , K and D satisfy the algebra 

[H X ,H 2 \ = [D,K} = 0, (4.9a) 
[K,Hi] =-[K,H 2 ] = -D, (4.9b) 
[D, Hi] = - [D, H 2 ] = -(H 1 + H 2 ), (4.9c) 

which differs from Eqs. ( |4.4|) by a change of sign in the right hand side. A possible realization 
of Eqs. ( |4.9D with a two component pseudopotential y = (yi,y 2 ) is 

Id Id 

H ^-^ + m- (410a) 
* = m - \h (410b) 

K = itf A (4.10c) 

D= Vh (410d) 



Hence, the pseudopotential equations ( |4.6|) take the form 

Vix = ^(ui + u 2 ) + -yl, (4.11a) 

y2x = ^(ui - u 2 ), (4.11b) 

Vit = \{.a\Uix + a 2 u 2x ) + -y 2 (l + otiui - a 2 u 2 ), (4.11c) 

o Z 

yit = ^(ociU lx - a 2 u 2x ). (4.11d) 

Equations ( |4. lla|) -( |4. llb|) produce 

ux = Ay lx + y 2x - 2y 2 2 , (4.12a) 

u 2 = Ay lx - y 2x - 2y 2 2 , (4.12b) 

which can be used to eliminate Ui and u 2 from Eqs. (|4.11c|) - (|4.11d| ). At this stage, let us seek 

travelling wave solutions of the quadratic model (|4.2[), i.e. solutions of the type U{ = Ui(£) 

where ^ = x + v t and i = 1,2. To this aim, let us start from Eq. (|4.11d ), which can be 
written as 

d ( 1 1 \ , . 

^7 y-vyi + -aiui - -a 2 u 2 j = 0. (4.13) 
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The integration of Eq. (|4 .13 ) gives 

1 



y\(0 



4(ai - o> 2 ) 



2c - 2vy 2 + 2(a 2 - u x )y 2 + (cm + a 2 )y' 2 



(4.14) 



with the help of Eqs. (|4.12|) , where Co is a constant of integration. Taking account of 
Eq. ([4.141) , from Eq. fl4.11c|) we get 



«i a 2 



y' 2 + v 



a 2 ai 



V 2 +2/2 



+ 



OiiOi 2 



a 2 ai 



(2c - 1; 



CpV 
OiiO! 2 



(4.15) 



where only the pseudopotential component y 2 is involved. This equation may be solved 
by the procedure exploited in Case III. Precisely, we could make a shift y — » y + y and 
fix the constants yo, cq in such a way that Eq. ( |4.15| ) takes the form Eq. ( |3.36| ). However, 
for simplicity, here we adopt a different approach: we take Co = and fix v by matching 
Eq. ( |3.36|) . Thus, we obtain 



aia 2 (a 2 - ai) 



\ (2a 2 — 3ai)(3a2 — 2ai) 



(4.16) 



which is real when 



Oil < gtt2, 



or 



Oi 2 < Oil < -a 2 . 



(4.17) 



Just to furnish an explicit numerical example, let us consider ai = 2, a 2 = 3/2. For- 
mula (|4.16| ) gives v = 5 and Eq. ( |4.15| ) becomes 



^ 4Q ^ 

As we have seen in Case III, a two-parameter solution of Eq. ( |4.18| ) is 
y 2 (0 = exp (If) P ^ exp + k u 0, k ) , 



(4.18) 



(4.19) 



with ko and fci arbitrary constants. In particular, taking ko = (see the discussion at the 
end of Case III) we obtain 

2 



2/2(0 



49 
20 



and, from Eqs. (|4.12|) and Eqs. Qj.l4j ), 



49 

Ml = ^(l + T) 2 (13 + 7T), 
49 

M2 = -(1 + T) 2 (8 + 7T), 
15 



T = tanh 



_7_ 
12 



(x — x + 5£) 



(4.20) 

(4.21a) 
(4.21b) 
(4.21c) 



with Xq arbitrary constant. 
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V. THE GIERER-MEINHARDT MODELS 



In this Section we apply the prolongation analysis to the RD models of the Gierer- 
Meinhardt (GM) type. In the biological context, of a special interest is the cubic GM 
system, which describes the interplay between an activator field u(x,t) and its substrate 
counterpart w(x, t). This model is defined by the two equations 

u t = au xx + R 1 , (5.1a) 

w t = (3w xx + R 2 , (5.1b) 

with 

Ri = e(u 2 w — u), (5.2a) 

R 2 = A(l -u 2 w). (5.2b) 

The two fields diffuse with diffusion coefficients which are generally different. The concen- 
tration of the field u decays according to the — eu term, but is enhanced by the substrate w 
via the production term eu 2 w. On the other hand, the substrate is injected in the domain 
of the reaction with a constant rate A while its depletion is controlled by the same nonlinear 
reaction term. 

The two opposite fixed points are 

u — 0, w = At + on, u> t — (3uj xx . 

The former describes a scenario in which the activator field has reached a balance between 
catalyzation by substrate and substrate depletion. In the latter case the activator field is 
absent and the substrate arises due to the injection constant term with inhomogeneities 
damped in time according to the heat equation. 

Moreover, according to the general analysis, it could be interesting also to investigate 
the model obtained by replacing in Eqs. ( |5.2|) the cubic interaction u 2 w with the quadratic 
term uw. We shall refer to this system as the quadratic Gierer-Meinhardt model. 

A. The cubic GM model 

Let us start from the prolongation equations 

y x = F(u,w,y), (5.4a) 
y t = G(u, u x , w, w x , y), (5.4b) 

whose compatibility condition provides 

aF u = G Ux , (5.5a) 

(3F w = G WxJ (5.5b) 

F U R\ + F W R 2 + [F, G\ = G u u x + G w w x . (5.5c) 

From these equations one finds that G u u = G WxWx = and, therefore, 
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G = A(u, w, y)u x + B(u, w, y)w x + C(u, w, y), 



(5.6) 



where A, B and C are vector fields of integration. Putting back this result into Eq. ( |5.5c|) 
and equating to zero the coefficients of the independent monomials 1, u x , w x , u x w x , u 2 x1 
we obtain 



A u = 0, 
B w = 0, 
A w + B u = 0, 
[F,A] =C U , 
[F, B] = C w , 
F U R X + F W R 2 + [F, C] = 0, 



(5.7a) 
(5.7b) 
(5.7c) 
(5.7d) 
(5.7e) 
(5.7f) 



from which 



A = a (y)w + ax(y), 
B = -a (y)u + a 2 (y), 

do, a\ and a2 being vector fields of integration. Moreover, Eqs. ( p.5a -5.5b) imply 



aF = G 

/3F W G Wx , 



Ma uw + a%u + h(w,y)), 



%(a u + K 



-aou + a 2 . 



Now, two cases have to be distinguished: if i) a + [3 ^ 0, then 

a = 0, 

a 

h = -a 2 w + a 3 , 

while if ii) a + (3 = 0, ao may be different from zero and 

h = —a 2 w + a 3 . 



(5.8a) 
(5.8b) 



(5.9) 



(5.10a) 
(5.10b) 



(5.11) 



First let us analyze the non-degenerate case i) which is the one relevant to proper RD 
systems [|TJ. 



1. The non- degenerate case: a + (3 ^ 
The equations to be satisfied are (see Eqs. 



[F,A]=C U , (5.12a) 
[F,B] =C W , (5.12b) 
+ F W R 2 + [F, C] = 0, (5.12c) 



where 
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A = 


a>i, 






F = 


a ( 






G = 





B = 02, (5.13a) 



a: 



aiu + —a 2 w + a 3 j , (5.13b) 
r + a 2 w x + C. (5.13c) 
Requiring that C uw = C wu (see Eqs. ( |5. 12"a] - |5. 12b|) ) , we find 

[oi,a 2 ]=0, (5.14) 

and 

1 1 

C = \au 03I it [02, atel w + 04- (5.15) 

a a 



Expanding Eq. ( 5.12c ) and collecting the coefficients of the monomials 1, u, w, u 2 , w 2 , uw 
, u 2 w we are led to the following incomplete algebra 

ci2 = ^—-cli, (5.16a) 
aX 

[ai, [ai,a 3 ]) = 0, (5.16b) 

[03,04] = — eax, (5.16c) 

[ai,a 4 ] = - [o 3 , [oi,a 3 ]] + eoi, (5.16d) 
a 

[a 1; a 4 ] = A [a 3 , [ai,o 3 ]] . (5.16e) 

This algebra corresponds to Eqs. ( |2.19| ). The linear dependence between ai and a 2 is due 
to the cubic terms appearing in the evolution equations ( |5.1|) . The compatibility condition 
Vtx — Uxt together with the information encoded into the incomplete algebra gives 

ai {first evolution eq.} ~ a 2 {second evolution eq.} . (5-17) 

On the other hand, since a\ and a 2 are linearly dependent, this equation is not equivalent to 
the pair of evolution equations ( |5.1| ). However, some non trivial structure remains because 
the cubic terms in Ri and R 2 are the same monomial in the fields. If this had not been the 
case, then a\ and a 2 would have been zero with a corresponding trivial algebraic structure. 
In order to close the above incomplete algebra, let us set 

[ai,a 3 ] = a B . (5.18) 

If we introduce the parameter 



9 — OL / \ 

e = ^— r , 5.19 



and rescale the generators as follows 
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dea 2 

°i = ~ <o ba ^i, (5.20a) 



a 3 = ~A 3 , (5.20b) 

a 4 = -r^-A 4 , (5.20c) 
p — a 

«s = - J ea \, A 5 , (5.20d) 



we have the commutation relations 

[A 1; A 3 ] =A 5 , [A 1 ,A A ]=A 1 , (5.21a) 
[A 1; A 5 ] =0, [A 3 , A 4 ] = Ax, (5.21b) 
[A 3 ,A 5 ] =A U [A 5 ,A A ] = A 5 , (5.21c) 

in terms of {Ai} i = 1, ... 5. Equations (|5.21|) define a closed Lie algebra as one can check 
by computing all the non trivial cases via the Jacobi identity. Moreover, if we put 

A 6 = A A - A 5 , (5.22) 

then Eqs. (|5.21|) become 

[A 1 ,A 5 ]=0, [A 3 ,A 6 ] = 0, (5.23a) 

[A 1 ,A 3 ]=A 5 , [A 1 ,A 6 ]=A l , (5.23b) 

[M,M =A U [X,A 6 ]=A B , (5.23c) 

where we recognize the Lie algebra sim(2) of the similitude group in the plane 0. The 
operators A\ and A 5 are the generators of two orthogonal translations, A 3 is the generator of 
rotations in the plane and, finally, A 6 is the generator of isotropic scalings. This geometrical 
interpretation gives also the basic vector field realization of the algebra, i.e. 

A = 7T, A 5 = — , (5.24a) 
ox ay 

d d 

A 3 =x—- y—, (5.24b) 
ay ox 

d d 

A 6 = x— + y—, (5.24c) 
ox oy 

It is rather intriguing the fact that the above rescaling requires a ^ (3 since this is also a 
necessary condition for Turing instability to occur in these kind of systems. The point a = f3 
is singular, in the sense that the incomplete algebra (|2.19f) collapses to 

a\ = 0, [03,04] = 0, (5.25) 

which implies the trivial result 

y x = -a 3 , (5.26a) 
a 

Vt = 04, (5.26b) 

[03, 04] = 0. (5.26c) 
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2. The degenerate case: a + (3 = 



From Eqs. ( |5.7| ) we deduce 



F = 


CLqUW + d\U — d2W + d 3 , 


(5.27a) 


G = 


Au x + Bw x + C, 


(5.27b) 


A = 


a(a,oW + ax), 


(5.27c) 


B = 


a(—a u + a 2 ), 


(5.27d) 


C = 


a {uw [a,i, a 2 ] — u[a\, a 3 ] — w [a 2 , 03] + 04} . 


(5.27e) 



We observe that in Eq. ( |5.12c ) the coefficient of the monomial u 3 w must vanish and this 



gives a = 0. Now, by performing the scaling e — > ae, A — > aA, the resulting prolongation 
algebra is 

(5.28a) 
(5.28b) 
(5.28c) 
(5.28d) 
(5.28e) 
(5.28f) 
(5.28g) 
(5.28h) 

Here, in theory, the generators a\ and a 2 are not constrained to be linearly dependent. 
However, complete linearization requires that specific algebraic constraints have to be satis- 
fied. For a generic model, some mechanism must ruin the equivalence between the evolution 
equations and the pseudopotential formulation. In this case, as we promptly show, a linear 
dependence arises among the generators. 
From Eq. (|5.28h| ) we find 



[a 2 , [a 2 ,ai]] 


= 0, 


[a 2 , [a 2 ,a 3 ]] 


= 0, 


[ai, [ai,a 3 ]] 


= 0, 


[a 3 , [ai,a 2 ]] 


= 0, 


[ai,a 4 ] 


= [a 3 , [ai,a 3 ]] + ea x 


[a 2 , o 4 ] 


= [a 3 , [03,02]] , 


[a 3 , a 4 ] 


= Aa 2 , 


[ai, [ai,a 2 ]] 


= — eax — Aa 2 . 



e [ ai , a 2 ] = [a 2 , [ai, [ai, a 2 ]]] = (5.29a) 
= [ai, [a 2 , [ai, a 2 ]]] - [[a x , a 2 ] , [a 2 , a x ]] = 0, (5.29b) 



hence, still from Eq. (|5.28h| ) 



a 2 = —~ai, (5.30) 
A 



and the incomplete algebra turns out to close on the algebra sim(2). 



B. The quadratic GM model 

In a similar way, also the quadratic Gierer-Meinhardt model turns out to be non lin- 
earizable. All the details of the analogous computations are contained in Appendix IB. 
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VI. SPECIAL SOLUTIONS OF THE NON-DEGENERATE CUBIC GM MODEL 



In this Section we focus on the cubic GM model in the non-degenerate case and discuss 
the consequences of the pseudopotential formulation, which can be written as 



u + \vj] ^A x y + ^-A 3 y, 
A / at 



y x = Fy = 

y t = Gy = - [u + ^-w I a^A x y + I u + I fiA 5 y + uA 4 y, 



i< \ _ r A _ , ( , A \ 



a A 



where 



/i = 



f3ea 



V 



6(3 



(3 — a ' 



a\ J 

2 = (3 -a 
ea 2 



(6.1a) 
(6.1b) 

(6.2) 



An explicit representation of the Lie algebra sim(2) is obtained from the geometric trans- 
formations of the similitude group and is given by 



/0 1\ 


\0 / 



/ 1 l\ 

-10 
\ 0/ 



(6.3) 



Aa 



/ 


-1 



















-l 


• 


■ -MS 





:) 


V 


















(6.4) 



Since the third row is null, the third component of the pseudopotential is constant, y% = c. 
The other equations turn out to be 



V2xx + 



Vl = -0l£V2x, 

1 uc / e 

TT^V2 = — [u + -w 
[at,) 1 a V A 

e(3 



Vn = A*c \u + — w J - uy 2 . 
The compatibility condition y 2 ,txx = U2,xxt gives the equation 

dt (u + = ad 2 r ( u + ^-w ) + e(l — it 



A 



Aa 



If we write u and iu in terms of the pseudopotential y = y 2 we find 

/?-al 



u = 1 + 



-(/3D X -D t )y, 



w = 



6(3 c 

\a [ j9-al . „ „ . 

■-T 1 + - {aD x -D t )y 

6(3 6a c 



(6.5a) 
(6.5b) 

(6.5c) 



(6.6) 

(6.7a) 
(6.7b) 



where 
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D * = d l + nw- ( 6 - 8a ) 

D t = d t + v. (6.8b) 

In order to determine y, we must put these expressions into one of the genuinely nonlinear 
evolution equations. When this is done, one obtains 

n <«W)v = jk [ + t i 1 + i w m ")' i 1 + kr n ^) } • < 6 - 9 > 

where fl(rj) = r)D x — D t and in particular VL(f3) = (331 — d t . 

In Appendix y we discuss the role of the Casimir operator C = A\ + A\ of the Euclidean 
subalgebra of sim(2). In particular, we show that if (u,w) is a solution of Eqs. ( |5.1|) and y 
a solution of Eqs. ( |6.1| ), then the new pseudopotential 

V^e-^Cy (6.10) 

is a new solution. 

Equation ( |6.9| ) will be applied below to find particular solutions (homogeneous and of 
the travelling wave type) to the cubic GM model. 

A. Homogeneous solutions 

Let us look for a class of solutions to Eqs. (|5.1|) assuming that 

u + i w = ( j ) (t), (6.11) 

A 

where <fi(t) denotes a given function of the time only. Then, from Eqs. ( |6.5|) we get 

u + 4 m = a? ^ + ^,_^... = 1 _^ + V^-l), (6.12) 



and therefore 



u = 1 0, (6.13a) 

w = -f-l+0+-0V (6.13b) 



e v e 

where <fi = dcf)/dt. The equations of motion (|5.1|) may be written as 



u t + \w t = e(l - u), (6.14a) 
A 

w t = \{l-u 2 w). (6.14b) 



The former is automatically satisfied by the Ansatz (|6.11 ), while the latter may be written 

as 
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0" = (1 - 00 (l + ^(1 - 00(1 - - 00) , (6-15) 

where 0' = d(f>/dr, with r = et. The (trivial) solutions corresponding to u = 1 and w = 
are respectively 

= 1 + |, = r-r o . (6.16) 

The problem to be solved is 

0" = (1 - 00 (1 + a(l - 00(1 - - 00) , « = -, (6-17) 



e 



with the initial conditions 



0'(O) = l-tt(O), (6.18a) 



0(0) = u(0) + T w(0). (6.18b) 
A 



Equation (|6.17|) can be elaborated by using the hodographic transformation 

1 



(6.19) 



r 1 + 0(0) 

In doing so, we obtain 

9' = -6 + 6 2 (a(f) - 2) + # 3 (a0 - 1 - a), (6.20) 

which gives for a = 1 

0' = 0(£0(0 + l)-l), (6.21) 
where 0' = d0/df , £ = - 2 and 

0'M = n I a L ov (6-22) 



It is interesting to consider the qualitative behavior of Eq. ( |6.21| ) with 6(0) = a > 0. This 
evolution problem 

0'(O = -0 + + £0 3 , (6.23a) 
6(0) = a > 0, (6.23b) 

has solutions which are decreasing until they eventually meet the nullcline T given by 6(6 + 
1)£ = 1. It seems reasonable to predict that all the solutions starting with a below some 
critical value a* decay exponentially at infinity whereas all the other solutions meet F and 
thereafter explode in a finite time. This scenario, described in Fig. El is confirmed by the 
exact solution of the evolution problem when only one of the nonlinear terms is present, 
namely 
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0' 

9' 



+ ^ = (l + 2^ + (l/^o - !)e 2€ 



-1/2 



(6.24a) 
(6.24b) 



A heuristic evaluation of a* is shown in Appendix |D]. Here we describe a notable approximate 
solution to the evoution problem ( |6.21| ). To this aim let us consider the equation 



o. 



0(0 + 1)" 2 (0+l) 

with £' = d^/dO. Now we remind the reader that the class of integrable equations 
allows the general integral 

/ = e ±? (e + /(0) + l). 
Actually, Eq. ( ggg ) can be associated which the evolution equation 

^ = Tj^j(* + /(*)). 



or 



d£ 
d£ 



+(£ + /), 



(6.25) 



(6.26) 



(6.27) 



(6.28) 



(6.29) 



which shows how integrability of Eq. ( |6.26| ) is just linearity in disguise. Anyhow, if we take 

1 



/(*) 



0(0 + 1) 



(6.30) 



and the minus sign in Eq. ( |6.28 ), then we find 



9' 



1) 



20 + 



_q W + i)_i] 



(6.31) 



This is a modified equation but with a modifying factor ^e+i^ which is bounded between 
1/2 and 1 when 9 > and which, therefore, may be expected to give minor changes in the 
solution. Another way of emphasizing the extra terms is that of writing 



d 



°4 



eme + i)-i]. 



(6.32) 



The integral of Eq. flPID is (see Eq. (|6727|) ) : 



In terms of 9 = 0(0) > we have 



0(0 + 1) 



+ 1 



(6.33) 



23 



1 

0o(0o+i; 



(6.34) 



and 



4 



l + £-/e«' 



(6.35) 



Some curves are displayed in Fig. [| The constant I is in the range (— oo, 1), all the solutions 
with I < are regular for £ > and decay exponentially at £ — > +oo. The particular value 
7 = corresponds to the separatrix between regular and singular solutions and starts at 
the golden ratio 9 = (\/5 — l)/2. We remark that Eq. ( |6.22p predicts singular solutions to 
be associated with periodic oscillating solutions approaching limiting closed curves in the 
(u, w) plane as t — > +oo. 

In the special case of the separatrix, all the remaining integrations may be performed. 
At I = we have 



0(0 



4 



1 + ^ 



9(<f>-2) 



+ 3 
- 1 : 



(6.36) 



and hence integrating Eq. ( |6.22| ) we find 
^0+^(0- 1)^0 + 30- 1 + log 
Moreover, from Eqs. ( |6.13| ) 



+ 1 + 



+ 3 



T - T . 



(6.37) 



U(T 



1 



^(0 + 3)7(0-1) 



+ 3)/(0-l) + l 



(6.38) 



that is 



- + u-l. 

u 



Finally, substitution from Eq. ( |6.39p into Eq. ( p.37| ) gives 

112 

- - H h log - = T - T . 

2 u u 



(6.39) 



(6.40) 



We remark that Eq. ( |6.40p can be exactly solved in terms of the Lambert W function [15 
which obeys the Schroder equation [JIB] 



W{x) + e w{x) = x. 



(6.41) 



Indeed, we have 



u 



W 



-e 2 



-t-tO 



(6.42) 
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The function w can be determined by Eqs. ( |6.13|) with the help of Eq. (|6.39| ), i.e. 



w = - - 1 = W ( le^ +T - T ° ) - 1. (6.43) 



u V2 

In Fig. [| we display the behavior of this special solution. To conclude this Section, we point 
out that the asymptotic behavior of W(x) for large x is W(x) ~ logx and therefore in this 
limit 

u(t) ~ (t-t + ~ -log 2^ , (6.44a) 

w(t) ~ r - r - ~ - log 2. (6.44b) 

We remark that this solution has been obtained from the approximate equation ( |6.31 ); how- 
ever, it solves asymptotically the exact equations ( |5.1|) because the extra terms in Eq. (|6.32|) 
vanish in the limit u — > (see Eqs. ( |6.13| ) and ( |6.19| )). 



B. Travelling waves 



In order to look for solutions of Eq. (|6.9| ) of the travelling wave type, we assume that the 
pseudopotential y is a function of the reduced variable £ — x + vt, where v is a constant. In 
terms of £ we have 

Q(a) = I3D 2 - vD, (6.45a) 
n((3) = aD 2 -vD-e, (6.45b) 

with D = This formalism suggests the particular speed 



v = ±PJ— ? , (6.46) 



as a critical value. Actually, in this case we have the factorization 



n {a) = a[ D±±—— lD T J—, (6.47a) 



with a common factor. Therefore, we can set 



and get a reduced equation in terms of the variable (p. Without loss of generality, let us 
focus on the particular set of constants 

a = 2, /3 = 1, e = 1. (6.49) 
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Then, v = ±1 and the relation between p and y becomes ip = (D^p l)y. Equation 
be written as 



(D±l/2) (D^l)D<p 



1 



1 - 2A (Dp - l) 2 ((£> ± 1/2) (p - 1) 



and the fields u, w are expressed in terms of p by 



u = 1 — D<^, 

io = -2A(l-(£>±l/2)y>) 



^9[) may 
(6.50) 



(6.51a) 
(6.51b) 



The fixed points of the cubic Gierer-Meinhardt model (see Eq. 
exact (trivial) solutions 



1 + 2A 



p = £ + const. 



IjD) correspond to the two 



(6.52) 



It is very interesting to look for a particular phenomenological meaning of the critical value 
expressed by Eq. ( |6.46| ). Indeed, our numerical simulations show that something special 
happens at that point (see Sec. [VI I| ). 

Apart from the critical speed, the general problem of finding a travelling wave solution to 
Eq. ( |6.9| ) may be studied by looking for interpolating solutions which at £ — > ±oo approach 
the two fixed points Q6.52 ). This is a difficult eigenvalue problem for a 4th order nonlinear 
ordinary differential equation. We shall exhibit semi-analytical approximations of these 



interpolating solutions in Sec. VII, devoted mainly to numerical results. 



VII. NUMERICAL SIMULATIONS 

In this Section we discretize the evolution equations of the cubic GM model and study 
them from a numerical point of view. This allows for a check of our analytical predictions 
and makes a step forward from our approximate solutions toward the unknown exact ones 
and their phenomenology. 

A. Homogeneous solutions 

Let us consider Eqs. Q5.1|) and ( |5.2|) in the case A = e (a = 1). By rescaling the time 
variable we can set 

A = e = l. (7.1) 

The couple of differential equations to be integrated numerically is 

u(t) = u 2 w — u, (7.2a) 
w(t) = 1 - u 2 w. (7.2b) 

The initial values w(0), w(0) are related by 

uj(0) = 2-u(0). (7.3) 
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Since the theoretical analysis carried out in Sec. [VI A| shows that oscillating trajectories 
emerge for 



u(Q) > 



a" 



1 + a* 



(7.4) 



where a* is a critical value determined in Appendix 0, we shall start with w(0) slightly above 
oi* /(l + a*). Using the following second order algorithm for the integration of equations 



in) M , °r/ (n) \ 



with 5 the discretization time step, we obtain the result pictured in Fig. 



(7.5a) 
(7.5b) 



B. Travelling waves 

The search for travelling wave solutions of the cubic GM model amounts to the solution 
of a 4th order nonlinear ordinary differential equation (a reduction of Eq. (|6.9|) ) with given 
boundary conditions corresponding to the two fixed points ( |5.3| ). From a biological point 
of view, these solutions interpolate between regimes dominated by the activator field or the 
substrate field and are characterized by a reaction front which separated the two regions. 
For a 2nd order equation, phase-space techniques permit geometrical proofs of the existence 
of such solutions. Here, dealing with a 4th order equation we shall build up approximate 
solutions and adopt them as particular initial conditions in the numerical integration. 

At the critical speed, the order is only 3 and we start from this case. The equation to be 
solved is Eq. ( |6.50| ). Without loss of generality, let us set A = 1. The two fixed points ( |6.52| ) 
are 

ip = £ + constant, <p = 3. (7.6) 

As we shall see, the above constant can be set to zero. Let us look for a perturbed solution 
around (p = £ by writing 

= £ + (7.7) 
and expanding Eq. ( ft.50 ) at the first order in 5. We find 

25"' - 5" - 5' = 0, (7.8) 



whose characteristic roots are 0, 1, —1/2. The same computation starting from 

<p = 3 + 5, (7.9) 

gives 

25"' -5" -5' -6 = 0, (7.10) 
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which has one characteristic root on the right complex half-plane and the other two on 
the left half-plane. A sensible perturbation must vanish at infinity. In order to have the 
greatest number of free constants we match the first fixed point at £ — > — oo and the other 
at £ — > +00. Hence, we construct our approximate solution as 

with 

<MO = £ + Ci + c 2 e«, (7.12a) 

(p + (g) = (c 3 sin v$ + c 4 cos «e~^ + 3, (7.12b) 

H = 0.366876, (7.12c) 

u = 0.520259. (7.12d) 

The constant C\ may be eliminated and set to zero by a translation in £. The free constants 
C2, C3 and C4 are determined (rather arbitrarily) by imposing the regularity condition 

v ,W (0) = v ,W (0)j n = 0) l,2. (7.13) 

The result is 

c 2 = 0.225361, (7.14a) 

c 3 = 0.398672, (7.14b) 

c 4 = -2.77464, (7.14c) 

which gives an approximate solution as regular as possible and asymptotically exact. 

This procedure may be extended to the more general case of an arbitrary speed. Since 
our interest is in the methodological viewpoint and in giving specific examples we keep the 
above parameters (a = 2, (5 = A = e = 1) but do not fix the speed. With these values, the 
travelling wave solutions of Eq. (|6.9|) obey 



2y vv - 'ivy" + (v 2 - l)y" + vy' = (7.15a) 
= -1 + (/ - vy' - l) 2 (2y" -vy' -y- 2), (7.15b) 



and the fields are expressed by 



u — 1 + vy' — y , (7.16a) 
w = -2-y -vy' + 2y". (7.16b) 

The two fixed points correspond to the solutions 

y=-3, y = -~. (7.17) 

v 

Perturbing around the first by setting y = —3 + 5 and expanding at the first order in 5 we 
obtain a linear differential equation with constant coefficients and characteristic polynomial 
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2p 4 - 3jA; + (V - iy + 1 = 0. (7.18) 

In Appendix [E] we show that for each value of the speed v the above equation has always 
two zeroes on the left complex half-plane and the other two on the right. Perturbing the 
second solution by setting y — —^/v + 5 we find 



2 v ± \Jv 2 + 8 
p(p — v)(2p — pv — 1) = 0, — >p = Q,v, , (7-19) 

and, for the same reasons as before, we are forced to match this boundary condition on the 
left at £ — ► — oo . Just to give a numerical result, the above procedure in the case v = 2 
provides the following approximate travelling wave 

with 

y _ = Cie ^ + C2e e(i+V3)/2_| (721a) 

y+ = -3 + (c 3 cos z/f + c 4 sin z/£) e ~^; (7.21b) 

where regularity up to the third order requires 

ci = 0.255476, (7.22a) 

c 2 = -0.690277, (7.22b) 

c 3 = -1.217923, (7.22c) 

c 4 = 2.565199. (7.22d) 

Turning to the numerical simulations, we fixed all the constants at the values a = 2.0, 
/3 = 1.0, e = A = 1.0. Assuming the equations of motion to be stabilized by the diffusion 
terms, we have discretized them by a first order Euler scheme. 

Concerning boundary conditions, they are easily dealt with at least in the case of trav- 
elling wave solutions. For generic positive speed v, the exact u wave has asymptotically 
constant values whereas the exact w wave tends to a constant on the right and to a linear 
function to the left. Using these informations we may impose the correct boundary con- 
ditions minimizing finite size effects. Anyway, we also checked independence of the results 
from the boundary conditions, at least when time is enough small for the size effects not to 
be relevant. For instance, if one exploits periodic boundary conditions and starts with the 
approximate (infinite volume) solution, then a perturbation is seen to arise at the boundary 
and the simulation must be stopped when it collides with the localized reaction front of the 
u wave. 

Now, let us examine the output of the numerical simulation. In Fig. |6| we display the 
evolution of the approximate wave with v = 2.0. As one can see, it rapidly settles down to 
a stable wave travelling with the desired speed. In Fig. |7] we repeat the simulation starting 
from the approximate wave with speed v = 1.0. In this case the approximate wave evolves 
toward an apparently stable travelling wave with internal oscillations. In Fig. |8|, where 
v = 0.5 is subcritical the same situation occurs. 
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These results suggest the following two scenarios below the critical wave speed: (i) there 
are no travelling waves or they are unstable; (ii) a stable travelling wave exists, but our 
approximate solutions are too crude and their evolution does not converge to it. Actually, 
in literature, we find examples of systems of this kind which exhibit the scenario (i) 0. 



VIII. CONCLUSIONS 

The results achieved in this paper show that the prolongation technique reveals a fruitful 
tool to investigate morphogenesis and autocatalysis modeling equations. 

We have dealt with a class of RD systems including reaction terms which are both 
quadratic and cubic in the fields. This class comprises well-known RD equations such as the 
Gierer-Meinhardt models. The prolongation approach allowed us to discover the existence 
of an algebraic structure inherent in this kind of models. It turns out to be the algebra 
associated with the similitude group. This feature could be important within a programme 
to settle up a systematic of RD models, which might be classified according to their algebraic 
properties. 

The comparison between special solutions drawn from theory and the corresponding 
ones obtained via numerical simulations is satisfactory. This represents an encouragement to 
extend the algebraic strategy to handle more complicated RD nonlinear evolution equations, 
keeping in mind cases in two-space and one-time dimensions. 

A direct search for particular solutions is possible as shown in [T7| where the socalled 



tanh method is utilized for a special cubic model. However, with this technique, calculations 
become rapidly cumbersome and are not suitable for model classification. Just to give an 
example, we considered the cubic model 

u t = au xx — u 2 w + Z\U + z 2 w + z 3 , (8.1a) 
w t = (3w xx — u 2 w + h\u + h 2 w + hz, (8.1b) 

and looked for a travelling wave with speed v. In the case a = f3 = v 2 z\ / (2z 2 ) , hi = z± + z 2 , 
h 2 = 0, h 3 = Zs/zx(zi + z 2 ), we were able to find the solution 

u = \fz\ + 2 2 tanh \-—yfz\ + z 2 ) , (8.2a) 
w = — + u, (8.2b) 



which, in the special case z 2 = 0, reproduces one of the solutions in |L7| . 

To conclude our comments, we notice that mathematical models are only a rough simpli- 
fication of a complex reality where the mechanisms of the biological, chemical and physical 
processes involved are often unknown. Notwithstanding, we think that an important role 
of these models is to get insight into possible interactions between specific processes and to 
suggest sometime new experiments |0 
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APPENDIX A: SOLUTION OF A SECOND ORDER DIFFERENTIAL 

EQUATION 

The solution of Eq. ( |3.36| ) can be achieved via a trick by Mittag-Lefner [0. In doing so, 
if Y(£) is a solution of the equation 

Y" = af - — a 2 Y + cY 2 , (Al) 
with a, c arbitrary constants, then 

H{0 = {Y'-\aY) 2 - 2 -cY\ (A2) 



fulfills the equation 



5 

Thus, we can find a constant Hq such that 

2 ~\ 2 2 ~ Q , /6 



H' = -aH. (A3) 



Y'-~aY) -~cY 3 = H exp^~ati). (A4) 
By making in ( |A4[ ) the change of variables 

F = exp(jU)^, -=^vf eXP (^)' (A5) 

we get 

tp' 2 {z) = 4:Lp 3 (z) — k, k constant, (A6) 

wich is solved by 

f(z)=V(z-z o ,0,k). (A7) 



APPENDIX B: PROLONGATION OF THE QUADRATIC MODEL 

In theory, quadratic reaction functions are compatible with linearizability. However, 
algebraic constraints have to be satisfied. As we shall see, the quadratic Gierer-Meinhardt 
model does not fall into the linearizable class, but showing this is not trivial. 
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a. The non- degenerate case: a + /3 ^ 

Repeating the same computation done in the case of the non-degenerate cubic model we 
find the following structure for the prolongation problem 

y x = F = -a x u + \a 2 w + a 3 , (Bla) 
a p 

y t = G = a x u x + a 2 w x - [a 1; a 3 ] u - [a 2 , a 3 ] w + a 4 , (Bib) 

The reaction equation 

F u Ri + F„,^ 2 + [F, - [at, a 3 ] u - [a 2 , a 3 ] w + a 4 ] = 0, (B2a) 
i^efW-tt), R 2 = \(l-uw), (B2b) 



leads to the following algebra 



[ai,a 2 ] 


= 0, 


[ai, [ai,a 3 ]] 


= 0, 


[a 2 , [a 2 ,a 3 ]] 


= 0, 


[ai,a 4 ] 


= eai + a [a 3 , [ai,a 3 ]] , 


[a 2 ,a 4 ] 


= p [a 3 , [a 2 ,a 3 ]] , 


[a 3 ,a 4 ] 


A 




- [ai, [02,^3]] 
a 


+ ^ [a 2 , [at, a 3 ]] = -ai 



A 

rf' 



(B3a) 
(B3b) 
(B3c) 
(B3d) 
(B3e) 

(B3f) 
(B3g) 



Let us scale the generators 

at/oc— > at, a 2 //3^a 2 , (B4) 

in order to cast the incomplete algebra into the simpler form 

[a u a 2 ] = 0, (B5a) 

[at, [at,a 3 ]] = 0, (B5b) 

[a 2 , [a 2 ,a 3 ]\ = 0, (B5c) 

[at, ^4] = ecu + oc [a 3 , [at, a 3 ]] , (B5d) 

[a 2 , a 4 ] = /3 [o 3 , [o 2 , a 3 ]] , (B5e) 

[a 3 , a 4 ] = -Aa 2 , (B5f) 

[ai, [a 2 , 03]] + [a 2 , [a u a 3 ]] = eai - Aa 2 . (B5g) 

Since [ai, a 2 ] = we can write 

[at, [a 2 , a*.]] = [a 2 , [ai, a fe ]] , VA;, (B6) 

and Eq. ( B5g|) may be rewritten as 
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[at, [a 2 , a 3 ]] = [a 2 , [at, a 3 ]] = -(eai - Aa 2 ). (B7) 



Moreover, from Eq. (P5d ) and Eq. ( |B5ej ) we find 



[a 2 , [at, o 4 ]] = a [a 2 , [a 3 , [ai, a 3 ]]] , (B8a) 
[ai, [a 2 , a 4 ]] = P [at, [a 3 , [a 2 , a 3 ]]] . (B8b) 



Then, using again Eq. ( p6|) we get 

a [a 2 , [a 3 , [at, a 3 )]] = (3[a x , [a 3 , [a 2 , a 3 ]]] , (B9) 

from which 

a [o 3 , [a 2 , [ai, a 3 ]]] - a [[ai, a 3 ] , [a 2 , a 3 ]] = (BlOa) 
= (3 [a 3 , [ai, [a 2 , a 3 ]]] - /3 [[a 2 , a 3 ] , [ai, a 3 ]] , (BlOb) 

or 

[[ai, a 3 ] , [a 2 , a 3 ]] = ^ - (e [01, a 3 ] - A [a 2 , a 3 ]), (Bll) 
zp + a 

by virtue of the Jacobi identity. Now if we substitute this result into the Jacobi identity 
among the three operators 

at, [01,03], [a 2 , a 3 ], (B12) 

we obtain [at, [a 2 , a 3 ]] = and consequently 

e ai = A a 2 . (B13) 

This linear dependence between at and a 2 constrains the incomplete algebra to close col- 
lapsing onto sim(2). In particular, complete linearizability is not achieved. 



b. The degenerate case: a + (3 = 
In the degenerate case, after proper rescaling, the incomplete algebra turns out to be 



[at, a 4 ] = [a 3 , [at, a 3 ]] + eat ~ Aao, (B14a) 

[at, [at, o 3 ]] = = [a 2 , [a 2 , a 3 ]] , (B14b) 

[a 2 ,a 4 ] = [a 3 , [a 3 ,a 2 ]] , (B14c) 

[a , 04] = -2 [a 3 , [ai, a 2 ]] + ea - eai - Aa 2 , (B14d) 

[ai, [ai,a 2 ]] = Aa , (B14e) 

[a 2 , [ai,a 2 ]] = ea , (B14f) 

[a 3 ,a 4 ] = Aa 2 , (B14g) 

[aojfli] =0, i = 1,2,3. (B14h) 
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It is easy to show that ao must vanish. In order to see this we write 

[ax, [a , a 4 ]] = -2 [ai, [a 3 , [ai, a 2 ]]] - A [ai, a 2 ] = 2 [[ai, a 2 ] , [ai, a 3 ]] - A [ai, a 2 ] , (B15) 
but the left hand side vanishes because 

[ax, [ao, 04]] = [a , [ai, a 4 ]] - [a 4 , [01, a ]] = [ao, [03, [ax, a 3 ]] + ea x - Aa ] = 0. (B16) 

Hence 

[[ai, a 2 ] , [ai, a 3 ]] = ^ [ai, a 2 ] , (B17) 



and 



2 

Aa = [ai, [ai, a 2 ]] = - [ax, [[ax, a 2 ) , [ax, o 3 ]j] = 
2 

= ^ {[[ai, a 2 ] , [ai, [a 1; a 3 ]]] - [[a 1; a 3 ] , [a x , [a u a 2 ]]]} = 
2 

= -- [[ai, a 3 ] , a A] = 0. (B18a) 

The algebra simplifies and becomes 

[ax, a 4 ] = [a 3 , [a x , a 3 ]] + eai, (B19a) 

[ai, [ai, a 3 ]] = = [a 2 , [a 2 , a 3 ]] , (B19b) 

[02,04] = [a 3 , [a 3 ,a 2 ]] , (B19c) 

[ai, [01,02]] =0, (B19d) 

[a 2 ,[a 1; a 2 ]] =0, (B19e) 

[a 3 , a 4 ] = Aa 2 , (B19f) 

[o 3 , [ai,a 2 ]] = --(eai + Aa 2 ). (B19g) 

Let us show that if the dimension of the space over which the above operator are defined is 
finite, then 

[a 1: [ax, a 2 ]] = = [a 2 , [a 1: a 2 ]] ^> [a x , a 2 ] = 0, (B20) 

from which, by using Eq. ( |B19gj ), we obtain the desired constraint 

eax + Xa 2 = 0. (B21) 

The proof of Eq. (|B20|) goes as follows: let A, B, X be a triple of linear operators acting on 
a finite dimensional vector space, satisfying 

[A, B] = X, [A, X] = 0, [B, X] = 0. (B22) 

Obviously TrX = 0. Furthermore, since 

Tr(A[B,C]) = Tr(B[C,A}), (B23) 
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B,X 



n-l 



0. 



we obtain, for n > 1, 

TrX n = TrpT 1 - 1 [A, B]) = Tr(A 
Thus, the characteristic polynomial of X is 

det(X - A) = (-A) dimX , 
and X has only the null eigenvalue. Let the canonical form of X be 

Mi \ 



V 



An J 



(B24) 



(B25) 



(B26) 



(where A n are n-dimensional Jordan matrices). The general form of the matrices which 
commute with (|B26| ) is block-rectangular where every block is a Toeplitz upper triangular 



matrix | 20fl . Then, it is easy to show that the condition [A, B] = X implies that all the A n 
blocks are 1-dimensional and therefore X = 0. 



APPENDIX C: THE SIMILITUDE GROUP AND CASIMIR OPERATORS 

Let us recall that given a Lie algebra C with generators T a 

[T a ,T p ] = clpT y , (CI) 
a Casimir operator C is a polynomial in {T a } which commutes with all the generators 



[C, T a ] = 0. 

If the Lie algebra is semisimple, the metric tensor 

9 a/3 = Ca\Cp T , 

is non singular and the possible Casimir operators are contained in the sequence 

T _ A „Pl rpa 1 rpa2 

1 2 — C Ql/3l C Q2/3 2 J 1 ' 

j _ J?2 J>3 3\ rpairpa 2 rpa 3 

J 3 — c ai/3i C a 2 /32 C a3/33 J 1 1 ' 



where 



(C2) 



(C3) 



(C4a) 
(C4b) 
(C4c) 



(C5a) 
(C5b) 



In the non semisimple case (for instance sim(2)), the independent Casimir operators must 
be built explicitly. Their number N is expressed by 



N = dim£ 



max rank||c£_aJ|, 

at,— ,a<m m c 



(C6) 
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according to the Beltrametti-Blasi theorem | 2l| . In the particular case of sim(2) we find 
N = and no Casimir operator exists. On the other hand, the Euclidean subalgebra 
generated by Ax, A 3 and A 5 admits a Casimir operator 



C = A\ + Al 

and since 

[C,A 4 ] = 2C, 

it is easy to show that given a solution y to the linearized problem 

y x = F x A x y + F 3 A 3 y, 

y t = G x A x y + G 5 A 5 y + uA^y, 



then the new pseudopotential 



-2vt 



(C7) 



(C8) 



(C9a) 
(C9b) 



(CIO) 



is another solution corresponding to the same (u, w) appearing in the functions F and G. 



APPENDIX D: HEURISTIC DETERMINATION OF a* 

Let us look for a solution to the equation 

6' = -6 + £6 2 + £6 3 , (Dl) 

with 

0(0) = a > 0, (D2) 

in the following form 

oo 

0(0 = E ^Pn(0, P(0 = a, P n <i(£) = 0, (D3) 

n=l 

where P n are functions to be determined. Choosing the integration constants in a sensible 
way, we obtain 



Pn(0 = e ( "- 1)€ ^ e-^H \ "jj P m (t)P n _ m (t) + "jj P^Pn-m-l \ ■ (D4) 



n— 1 n— 1 



m=l m=l,(=l 



The functions P n are polynomials, the first two being Pi = a, Pj = — a 2 (£ + 1) • The initial 
value 6(0) is an unknown function of a admitting the series expansion 

3 17 60^ 

0(0) = E ^(0) = a - a 2 + -a 3 - -a 4 - —a 5 + • • • = E (D5) 
This series shows a poor convergence near a ~ 1, but in terms of the mapped variable 
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we find 

f(0)=^-^-^ 4 + --- = E^ B - (D7) 
If we assume this series to have an infinite convergence radius and plot the function 




we see that 0(0) appears to be upper bounded by a finite constant a* whose numerical value 
can be estimated by truncating the infinite coefficients b' n . Taking the first 20 coefficients, 
we find 

a* = 0.756561. (D9) 



APPENDIX E: ANALYSIS OF THE ROOTS OF A QUARTIC EQUATION 

If Q(x) is a polynomial of degree N, then the number Z of roots on the left side of the 
imaginary axis is given by 



2 2tt J-oo Q{z) 



dy. (El) 



z=iy 



In the case Q(x) = 2x 4 — 3vx 3 + {v 2 — l)x 2 + 1 we have 



Z = 2 + J- r H(y)dy, (E2) 

where 

H(y) = 4- I Arctan— 1 . (E3) 

Kyj dy \ 2y 4 + (1 - v 2 )y 2 + 1 J K ' 

Since Arctan(- ■ ■) — > as |y| — > we can conclude that Z = 2. The Arctan argument is odd 
and its possible singularities give contributions which cancel in pairs. 
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FIGURES 



FIG. 1. Exact solutions for the u, w fields in the four cases corresponding to Eqs. ( 3.23| ), ( [3.30D 



( 3.38| ) and ( 3.4CQ . The numerical values of the parameters are respectively i) v = 1 , ax = 1 , a 2 = 



, Ai = 1 , A 2 = -1 , iii = 1 , C = 1; ii) v = 1 > «i = 1 i «2 = -1 , Ai = 1 , A 2 = 1 , \i\ = 1 , C = i; 
iii) v = 1 , ai = 2 , a; 2 = 1 , Ai = 1 , A 2 = —2 , /zi = 1 , £ = 0; iv) *y = 1 , a± = 1 , a.2 = — 1 , 
Ai = 1 , A 2 = 1 , Hi = 1 , C = !, Pi = !• 

FIG. 2. Qualitative behaviour of Eq. ( |6.21| ). 
FIG. 3. Analytical integrals of Eq. (6.31). 



FIG. 4. Solution corresponding to Eqs ( 6.42; ) and (|6.43| ). 



FIG. 5. Homogeneous solutions obtained from numerical integration of Eqs. (ffj^). We remark 
that from the value of a* given in Eq. ( |D9| ) we get a*/(l + a*) ~ 0.431. 



FIG. 6. Evolution of the approximate travelling wave with v = 2.0. The values of the other 
parameters are a = 2, j3 = 1, e = 1, A = 1, Ax = 0.05, At = 0.00025, Ax and At being the space 
and time discretization steps. 

FIG. 7. Like Fig. §, but v = 1.0. 
FIG. 8. Like Fie. I, but v = 0.5. 
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